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DEVELOPMENT AND APPLICATION OF BENCHMARK EXAMPLES FOR MODE II 
STATIC DELAMINATION PROPAGATION AND FATIGUE GROWTH PREDICTIONS 

Ronald Krueger* 

ABSTRACT 

The development of benchmark examples for static delamination propagation and 
cyclic delamination onset and growth prediction is presented and demonstrated for a 
commercial code. The example is based on a finite element model of an End-Notched 
Flexure (ENF) specimen. The example is independent of the analysis software used and 
allows the assessment of the automated delamination propagation, onset and growth 
prediction capabilities in commercial finite element codes based on the virtual crack closure 
technique (VCCT). First, static benchmark examples were created for the specimen. 

Second, based on the static results, benchmark examples for cyclic delamination growth 
were created. Third, the load-displacement relationship from a propagation analysis and the 
benchmark results were compared, and good agreement could be achieved by selecting the 
appropriate input parameters. Fourth, starting from an initially straight front, the 
delamination was allowed to grow under cyclic loading. The number of cycles to 
delamination onset and the number of cycles during delamination growth for each growth 
increment were obtained from the automated analysis and compared to the benchmark 
examples. Again, good agreement between the results obtained from the growth analysis 
and the benchmark results could be achieved by selecting the appropriate input parameters. 

The benchmarking procedure proved valuable by highlighting the issues associated with 
choosing the input parameters of the particular implementation. Selecting the appropriate 
input parameters, however, was not straightforward and often required an iterative 
procedure. Overall the results are encouraging, but further assessment for mixed-mode 
delamination is required. 


1. INTRODUCTION 

Over the past two decades, the use of fracture mechanics has become common practice to 
characterize the onset and growth of delaminations. In order to predict delamination onset or 
growth, the calculated strain energy release rate components are compared to interlaminar fracture 
toughness properties measured over a range from pure mode I loading to pure mode II loading. 

The virtual crack closure technique (VCCT) is widely used for computing energy release rates 
based on results from continuum (2D) and solid (3D) finite element (FE) analyses and to supply the 
mode separation required when using the mixed-mode fracture criterion [1, 2]. The virtual crack 
closure technique was recently implemented into several commercial finite element codes. As new 
methods for analyzing composite delamination are incorporated into finite element codes, the need 
for comparison and benchmarking becomes important since each code requires specific input 
parameters unique to its implementation. These parameters are unique to the numerical approach 
chosen and do not reflect real physical differences in delamination behavior. 

An approach for assessing the mode I, and mixed-mode I and II, delamination propagation 
capabilities in commercial finite element codes under static loading was recently presented and 
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demonstrated for VCCT for ABAQUS® 1 [3] as well as MD Nastran™ and Marc™ 2 [4], First, 
benchmark results were created manually for full three-dimensional finite element models of the 
Double Cantilever Beam (DCB) and the Single Leg Bending (SLB) specimen. Second, starting 
from an initially straight front, the delamination was allowed to propagate using the automated 
procedure implemented in the finite element software. The approach was then extended to allow the 
assessment of the delamination fatigue growth prediction capabilities in commercial finite element 
codes [5]. As for the static case, benchmark results were created manually first. Second, the 
delamination was allowed to grow under cyclic loading in a finite element model of a commercial 
code. In general, good agreement between the results obtained from the propagation and growth 
analysis and the benchmark results could be achieved by selecting the appropriate input 
parameters. Overall, the results were encouraging but showed that additional assessment for mode 
II and other mixed-mode delamination cases are required. 

The objective of the present study was to create additional benchmark examples, independent of 
the analysis software used, which allows the assessment of the static delamination propagation as 
well as the onset and growth prediction capabilities in commercial finite element codes. For the 
simulation of mode II fracture, the three-point bending End Notched Flexure (ENF) specimen was 
selected as shown in Figure 1. Dimensions, layup and material properties were taken directly from a 
related experimental study [6]. 

Static benchmark results were created based on the approach developed earlier [3], using two- 
dimensional finite element models for simulating the ENF specimens with different delamination 
lengths ciq. For each delamination length modeled, the load, Q, and the displacement, w, (center 
deflection) were monitored. The mode II strain energy release rate, Gu, was calculated for a fixed 
applied load. It is assumed that the delamination propagates when computed energy release rate, Gu, 
reaches the fracture toughness Gu c . Thus, critical loads and critical displacements for delamination 
propagation were calculated for each delamination length modeled. From these critical 
load/displacement results, benchmark solutions were created. It is assumed that the 
load/displacement relationship computed during automatic propagation should closely match the 
benchmark cases. 

Benchmark cases to assess the growth prediction capabilities were created based on the finite 
element models of the ENF specimen used for the static benchmark case. First, the number of cycles 
to delamination onset, No, was calculated from the mode II fatigue delamination growth onset data 
of the material [6]. Second, the number of cycles during delamination growth, ANg, was obtained 
incrementally from the material data for mode II fatigue delamination propagation [6] by using 
growth increments of Aa= 0.1 mm. Third, the total number of growth cycles, Ng, was calculated by 
summing over the increments ANg ■ Fourth, the corresponding delamination length, a, was 
calculated by summing over the growth increments Aa. Finally, for the benchmark case where 
results for delamination onset and growth were combined, the delamination length, a, was 
calculated and plotted versus an increasing total number of load cycles Nj=Nd+Ng- 

After creating the benchmark cases, the approach was demonstrated for the commercial finite 
element code ABAQUS®. Starting from an initially straight front, the delamination was allowed to 
propagate under static loading or grow under cyclic loading based on the algorithms implemented 
into the software. Input control parameters were varied to study the effect on the computed 
delamination propagation and growth. The benchmark enabled the selection of the appropriate 


1 ABAQUS® is a product of Dassault Systemes Simulia Corp. (DSS), Providence, RI, USA 

2 MD Nastran™ and Marc™ are manufactured by MSC. Software Corp., Santa Ana, CA, USA. NASTRAN® is a 
registered trademark of NASA. 
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input parameters that yielded good agreement between the results obtained from the growth 
analysis and the benchmark results. Once the parameters have been identified, they may then be 
used with confidence to model delamination growth for more complex configurations. 

In this paper, the development of the benchmark cases for the assessment of the static 
delamination propagation as well as the onset and growth prediction capabilities are presented. 
Examples of automated propagation and growth analyses are shown, and the selection of the 
required code specific input parameters are discussed. 


2. METHODOLOGY 

For the current numerical investigation, the three-point End-Notched Flexure (ENF) specimen, 
as shown in Figure 1, was chosen since it is simple, only exhibits the mode II opening fracture mode 
and had been used previously to develop an approach to assess the quasi-static delamination 
propagation simulation capabilities in commercial finite element codes [4], The methodology for 
delamination propagation, onset and growth was applied to the ENF specimen to create the 
benchmark example [7, 8 ]. For the current study, an ENF specimen made of IM7/8552 
graphite/epoxy with an unidirectional layup, [0]24, was modeled. The material, layup, overall 
specimen dimensions including initial crack length, a , were identical to specimens used in a related 
experimental study [ 6 ]. The material properties are given in Tables I and II. 

2.1 Static fracture toughness 

A quasi-static mixed-mode fracture criterion is determined first, since the parameters are 
required input to the ABAQUS® VCCT solution as discussed in the appendix. The mixed-mode 
fracture criterion for a material is determined by plotting the interlaminar fracture toughness, G c , 
versus the mixed-mode ratio, Gu/Gt as shown in Figure 2. The fracture criterion is generated 
experimentally using pure Mode I (Gn/Gj=0) Double Cantilever Beam (DCB) tests, pure Mode II 
(Gn/Gj=l) End-Notched Flexure (ENF) test (as shown in Figure 1), and Mixed Mode Bending 
(MMB) tests of varying ratios of G/ and Gy/. For the material used in this study, the experimental 
data (open symbols) and mean values (filled symbols) are shown in Figure 2. The different symbols 
(circles and squares) denote different data sources for the material [ 6 , 9]. A 2D fracture criterion 
was suggested by Benzeggah and Kenane [10] using a simple mathematical relationship between G c 
and Gji/Gt 


G r -G, I + (G„ c -Gj' 


( n V 

G ii 

k G t ) 


( 1 ) 


In this expression, G/ e and Gu c are the experimentally-determined fracture toughness data for mode I 
and II as shown in Figure 2. The factor 17 was determined by a curve fit using the Levenberg- 
Marquardt algorithm in the KaleidaGraph™ graphing and data analysis software [11]. 

2.2 Fatigue delamination growth onset 

The number of cycles to delamination onset, No, can be obtained from the delamination onset 
curve plotted in Figure 3 [ 6 , 12], The onset curve (solid green line) is a power law fit 

G = m 0 ■ Np‘ (2) 
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of the experimental data (open, green circles) obtained from an ENF test using a draft standard for 
delamination growth onset [6]. 


2.3 Fatigue delamination growth 

The number of cycles during delamination growth, No, can be obtained from the fatigue 
delamination propagation relationship (Paris Law) plotted in Figure 4 [6, 12]. The delamination 
growth rate (solid black line) can be expressed as a power law function 


da 

dN 


G n 


( 3 ) 


where da/dN is the increase in delamination length per cycle and G max is the maximum energy 
release rate at the front at peak loading. The factor c and exponent n were obtained by fitting the 
curve to the experimental data obtained during the ENF tests [6]. The critical energy release rate or 
fracture toughness, Gn c , was included in the plot of Figure 4 (blue solid vertical line). Since 
composites do not exhibit the same threshold behavior commonly observed in metals, a cutoff 
value, G t h, was chosen below which delamination growth was assumed to stop (green solid vertical 
line) [6], 


3. FINITE ELEMENT MODELING 
3.1 Model description 

An example of a two-dimensional finite element model of an End-Notched Flexure (ENF) 
specimen with boundary conditions is shown in Figure 5. The specimen was modeled with solid 
plane strain elements (CPE4, CPE4I) and solid plane stress elements (CPS4) in ABAQUS® 
Standard 6.9EF and 6.10. Along the length, all models were divided into different sections with 
different mesh refinement as shown in Figure 5a. The ENF specimen was modeled with six 
elements through the specimen thickness ( 2h ) as shown in the detail of Figure 5b. The resulting 
element length at the delamination tip was Aa=1.0 mm. Finer meshes, resulting in Aa=0.5 mm and 
Aa=0.25 mm were also generated as shown in Figures 5c and d, respectively. 

The plane of delamination was modeled as a discrete discontinuity in the center of the specimen. 
For the analysis with ABAQUS® 6.9EF and 6.10, the models were created as separate meshes for 
the upper and lower part of the specimens with identical nodal point coordinates in the plane of 
delamination [13]. Two surfaces (top and bottom surface) were defined to identify the contact area 
in the plane of delamination as shown in Figure 5b. Additionally, a node set was created to define 
the intact (bonded nodes) region. 

Examples of three-dimensional finite element models of the ENF specimen are shown in 
Figures 6 and 7. Along the length, all models were divided into different sections with different 
mesh refinement. A refined mesh was used in the center of the ENF specimen as shown in the detail 
of Figure 6b. Across the width, a uniform mesh (25 elements) was used to avoid potential problems 
at the transition between a coarse and finer mesh [3-5]. Through the specimen thickness ( 2h ), six 
elements were used as shown in the detail of Figure 6b. The resulting element length at the 
delamination tip was Aa=1.0 mm. The specimen was modeled with solid brick elements (C3D8I), 
which had yielded excellent results in a previous studies [3.4], A refined mesh with an increased 
number of elements in the length direction, resulting in Aa=0.5 mm, was also generated (not 


4 



shown). For an additional model, the width direction was also refined (50 elements) resulting in the 
model shown in Figure 7a. 

An example of a continuum shell element (SC8R) finite element model of the ENF specimen is 
shown in Figure 7b. Typically, the continuum shell elements in ABAQUS® are used to model an 
entire three-dimensional body. Unlike conventional shells, which model a reference surface, the 
SC8R elements have displacement degrees of freedom only, use linear interpolation, and allow 
finite membrane deformation and large rotations and, therefore, are suitable for geometric nonlinear 
analysis. The continuum shell elements are based on first-order layer-wise composite theory and 
include the effects of transverse shear deformation and thickness change [13]. In the x-y plane, the 
models had the same fidelity as the models made of solid brick elements C3D8I shown in Figures 6 
and 7a, resulting in an element length at the delamination tip Aa= 1.0 mm, and Aa=0.5 mm, 
respectively. In the z-direction, only one element was used to model the thickness of each specimen 
arm. These less-refined models were used to study the effect on performance (CPU time), computed 
load/displacement behavior and growth prediction in comparison with the more refined solid brick 
(C3D8I) models shown in Figures 6 and 7a. 

3.2 Static delamination propagation analysis 

For the automated delamination propagation analysis, the VCCT implementation in ABAQUS® 
Standard 6.9EF and 6.10 were used. The plane of delamination in three-dimensional analyses is 
modeled using the existing ABAQUS®/Standard crack propagation capability based on the contact 
pair capability [13]. Additional element definitions are not required, and the underlying finite 
element mesh and model does not have to be modified [13]. The implementation offers a crack and 
delamination propagation capability in ABAQUS®. It is implied that the energy release rate at the 
crack tip is calculated at the end of a converged increment. Once the energy release rate exceeds the 
critical strain energy release rate (including the user-specified mixed-mode criteria as shown in 
Figure 2), the node at the crack tip is released in the following increment, which allows the crack to 
propagate. To avoid sudden loss of stability when the crack tip is propagated, the force at the crack 
tip before advance is released gradually during succeeding increments in such a way that the force is 
brought to zero no later than the time at which the next node along the crack path begins to open 

[13]. 

In addition to the mixed-mode fracture criterion, VCCT for ABAQUS® requires additional input 
for the propagation analysis. If a user specified release tolerance is exceeded in an increment 
0 G-G c )/G c > release tolerance, a cutback operation is performed which reduces the time 
increment. In the new smaller increment, the strain energy release rates are recalculated and 
compared to the user specified release tolerance. The cutback reduces the degree of overshoot and 
improves the accuracy of the local solution [13]. A release tolerance of 0.2 is suggested in the 
handbook [13]. To help overcome convergence issues during the propagation analysis, ABAQUS® 
provides: 

• contact stabilization, which is applied across only selected contact pairs and used to 
control the motion of two contact pairs while they approach each other in multi- 
body contact. The damping is applied when bonded contact pairs debond and move 
away from each other [13] 

• automatic or static stabilization which is applied to the motion of the entire model 
and is commonly used in models that exhibit statically unstable behavior such as 
buckling [13] 

• viscous regularization, which is applied only to nodes on contact pairs that have just 
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debonded. The viscous regularization damping causes the tangent stiffness matrix 
of the softening material to be positive for sufficiently small time increments. 
Viscous regularization damping in VCCT for ABAQUS® is similar to the viscous 
regularization damping provided for cohesive elements and the concrete material 
model in ABAQUS®/Standard [13]. Further details about the required input parameters 
are discussed in the appendix. 

For automated propagation analysis, it was assumed that the computed behavior should closely 
match the benchmark results created below. For all analyses, the elastic constants (given in Table I) 
and the input to define the fracture criterion (given in Table II) were kept constant. The following 
items were varied to study the effect on the automated delamination propagation behavior during 
the analysis: 

• The release tolerance ( reltol) was varied. 

• The input for contact stabilization (cs) was varied. 

• The input for viscous regularization ( damv ) was varied. 

• Two- and three-dimensional models with different types of elements were used. 

• Models with different element length at the crack tip/delamination front, Aa, were used. 

3.3 Fatigue delamination onset and growth analysis 

For the automated delamination onset and growth analysis, the low-cycle fatigue analysis in 
ABAQUS® Standard 6.9EF and 6.10 was used to model delamination growth at the interfaces in 
laminated composites [13, 14]. A direct cyclic approach is part of the implementation and provides 
a computationally effective modeling technique to obtain the stabilized response of a structure 
subjected to constant amplitude cyclic loading. The theory and algorithm to obtain a stabilized 
response using the direct cyclic approach are described in detail in reference 14. Delamination onset 
and growth predictions are based on the calculation of the strain energy release rate at the 
delamination front using VCCT. To determine propagation, computed energy release rates are 
compared to the input data for onset and growth from experiments as discussed in the methodology 
section. The implementation is set up to release at least one element at the interface after the loading 
cycle is stabilized [13, 14]. 

For automated delamination onset and growth analyses, it was assumed that the computed 
behavior should closely match the benchmark results created below. For all analyses, the elastic 
constants (given in Table I), the input to define the fracture criterion (given in Table II), and the 
parameters for delamination onset and delamination growth (Paris Law) (given in Table II) were 
kept constant. The parameters to define the load frequency (/= 5 Hz), the load ratio (/?=0. l ) as well 
as the minimum and maximum applied displacement ( w min and w max ) were also kept constant 
during all analyses. A Fourier series is used during the execution of the ABAQUS® Standard to 
approximate the periodic cyclic loading. Based on previous results, it was decided to use 50 
Fourier terms to approximate the periodic cyclic loading [5]. Note, that number of Fourier terms 
is required input in addition and independent of the definition of the applied cyclic loading. Further 
details about the required input parameters are discussed in the appendix where a sample input file 
is also provided. The following items were varied to study the effect on the predicted delamination 
onset and growth behavior during the analysis: 

• The input to define the cyclic loading was studied. The starting time, to, which causes a 
phase-shift was varied. 

• Two- and three-dimensional models with different types of elements were used. 
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For three-dimensional models, the number of terms used to define a Fourier series was 
reduced to decrease computation time. 

Models with different element length at the crack tip/delamination front, Aa, were used. 


4. DEVELOPMENT OF THE STATIC BENCHMARK CASE 

The static benchmark case was created based on the approach developed earlier [3]. Two- 
dimensional finite element models simulating ENF specimens with 15 different delamination 
lengths cio were created (25.4 mm< ao<162 mm). For each delamination length modeled, the load, 
Q, and displacement, w, were monitored as shown in Figure 8 (colored lines). Using VCCT, the 
mixed-mode strain energy release rate components were computed for the applied load 0=2000 N 
as shown in Figure 8. As expected, the results were predominantly mode II. Therefore, a failure 
index Gn/Gn c was calculated by correlating the results with the mode II fracture toughness, Gu c , of 
the graphite/epoxy material. It is assumed that the delamination propagates when the failure index 
reaches unity. Therefore, the critical load, Q crit , can be calculated based on the relationship between 
load, Q, and the energy release rate, G [15], 


Q 2 , dCp 
2 dA 


(4) 


In equation (4), Cp is the compliance of the specimen, and dA is the increase in surface area 
corresponding to an incremental increase in load or displacement at fracture. The critical load, Q a -it, 
and critical displacement, w crit , were calculated for each delamination length modeled 


G„ = Q 

Guc Qcrit 



(5) 


and the results were included in the load/displacement plots as shown in Figure 9 (solid red circles). 
These critical load/displacement results indicated that, with increasing delamination length, less 
load is required to extend the delamination. For the first ten delamination lengths, ao, investigated, 
the values of the critical displacements also decreased at the same time. This means that the ENF 
specimen exhibits unstable delamination propagation under load as well as displacement control in 
this region. The remaining critical load/displacement results pointed to stable propagation. 

From these critical load/displacement results (solid red line), two benchmark solutions can be 
created as shown in Figure 10. During the analysis, either prescribed displacements, w, or nodal 
point loads, Q, are applied. For the case of prescribed displacements, w, (dashed blue line), the 
applied displacement must be held constant over several increments once the critical point (Qa-it, 
v C rit ) is reached, and the delamination front is advanced during these increments. Once the critical 
path (solid black line) is reached, the applied displacement is increased again incrementally. For the 
case of applied nodal point loads (dashed red line), the applied load must be held constant while the 
delamination front is advanced during these increments. Once the critical path (solid black line) is 
reached, the applied load is increased again incrementally. 

The benchmark result for prescribed displacements may also be visualized by plotting the 
prescribed displacements, w, at delamination growth onset versus the increase in delamination 
length, a*, as illustrated in Figure 11. For the case of applied nodal point loads, the benchmark 
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result may be visualized by plotting the applied loads, Q, at delamination growth onset versus the 
increase in delamination length, a* as illustrated in Figure 12. It is assumed that the 
load/displacement, load/delamination-length or applied displacement/delamination-length 
relationship computed during automatic propagation should closely match the benchmark cases. 


5. CREATING A BENCHMARK EXAMPLE FOR GROWTH PREDICTION 

For the development of the benchmark case for delamination onset from an initial flaw and 
subsequent growth under cyclic loading, guidance was taken from test results for mode II 
interlaminar fracture toughness and fatigue characterization [6]. In the report, a series of tests are 
documented which were conducted under load control with a maximum load, Q max , which caused 
the energy release rate at the front, Gu mca , to reach values of 50% of Gn c 

Gllma*_ = Q 5 (6) 


The maximum load, Q 50 , max-, and the maximum displacement, wso.max, were calculated using the 
known quadratic relationship between energy release rate and applied load or displacement 


Umax 


' lie 


Ql 


Q, 


2 

crit 



( 7 ) 


Qmax = Qcrit ^ , W ^ = W crU VOA 


( 8 ) 


where Q cn , and w a -u are the critical values. For the current study, a critical energy release rate 
Gnc = 0.78 kJ/m was used, and the critical values, Q crit , (solid red horizontal line) and, w crih (solid 
red vertical line) were obtained from the benchmark for static delamination propagation (discussed 
above) shown in the load-displacement plot in Figure 13. The calculated maximum load, Qso.max, 
and calculated maximum displacement, w 50 , max, are shown in Figure 13 (dashed red line) in 
relationship to the static benchmark case (solid grey circles and dashed grey line) mentioned above. 
During constant amplitude cyclic loading of an ENF specimen under load control, the applied 
maximum load, (9jft OT < K =1079 N, is kept constant while the displacement increases with increasing 
delamination length (horizontal dashed red line). For simulations performed under displacement 
control, the applied maximum displacement, wjo jOTax =1.0 mm, is kept constant while the load 
decreases as the delamination length increases (vertical dashed red line). In the report on mode II 
interlaminar fracture toughness and fatigue characterization [6], tests are also documented that were 
performed at 40%, 30% and 20% of Gnc The corresponding maximum loads and the maximum 
displacements were calculated using equation (8) and were included in the plot of Figure 13 
(Qffl.maxi Q 30 ,max> Q 20 ,max ~ horizontal lines, W 40 : maxi w 30 , max ? W 20 , max ~ vertical lines). 

For each of the 15 finite element models representing 15 delamination lengths (25.4 mm< 
ao<162 mm) as shown in Figure 8, the energy release rate corresponding to an applied load, 
Q 50 ,max = 1079 N, was calculated using equation (7). The energy release rate first increased with an 
increase in delamination length, a* as shown in Figure 14 (open red circles and dashed red line). 
After reaching a peak, the energy release rate decreased with increasing delamination length. 
Delamination growth was assumed to become unstable once the calculated energy release rate 
passed the fracture toughness value Gn c (solid blue line in Figure 14). Later, delamination growth 
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was assumed to become stable again after the calculated energy release rate dropped below the 
fracture toughness value Gn c . Additionally, the energy release rate dependence on the crack length 
was calculated for Qw.max, Qso.max, and Q2o,max, and the results were included in the plot of Figure 14 
(dashed lines with open symbols). The static benchmark case (solid grey circles and dashed grey 
line in Figure 14), where the delamination propagates at constant Gn c (solid blue line in Figure 14) 
was included for comparison. Also included was the cutoff value, G t h, (green solid horizontal line). 

For each of the 15 finite element models representing 15 delamination lengths (25.4 mm< 
a o<162 mm) as shown in Figure 8, the energy release rate corresponding to an applied 
displacement w 50imax = 1 .0 mm was also calculated using equation (7). The energy release rate first 
increased with an increase in delamination length, a*, as shown in Figure 15 (solid red circles and 
dashed red line). After reaching a peak, the energy release rate decreased with increasing 
delamination length. Delamination growth was assumed to stop once the calculated energy release 
rate dropped below the cutoff value, G t h, (green solid horizontal line). Additionally, the energy 
release rate dependence on the crack length was calculated for W 4 o, ma x, wso.max, and w 20 , max, and the 
results were included in the plot of Figure 15 (dashed lines with solid symbols). The static 
benchmark case (solid grey circles and dashed grey line in Figure 15), where the delamination 
propagates at constant Gn c (solid blue line in Figure 15) was included for comparison. 

The ENF fatigue characterization tests were performed at a load ratio R= 0.1 [6]. The 
corresponding minimum load, Q min , and minimum displacement, w min „ were calculated 

R= 0mun_ = o.l => Q = 0.1 Q and w =0.1-w (9) 

^ x^min xZ-max min max V / 

O w 

*^max max 

Further, the tests were performed at a frequency f=5 Hz [6]. A graphical representation of the cyclic 
fatigue loading is plotted in Figure 16 for the example of displacement control Wio imax =\ .0 mm. The 
applied displacement, w, is represented as a function of time, t 

w = [A 0 + 5i-sinm(t-t 0 )]-w m(K (10) 

where w mca is the maximum displacement. The constants Aq= 0.55, Bi= 0.45, the circular frequency 
(o=10k= 3I.42 and the starting time fo=0.05 are calculated from load ratio R= 0.1 and the frequency 
f= 5 Hz for testing. The resulting equation to calculate the applied displacement, w, is shown in 
Figure 16. 

5.1 Fatigue delamination growth onset 

The number of cycles to delamination onset, No, may be obtained by solving equation (2) for 

N d . 


( ^ X, 


G = m 0 ■ Np 1 


N d = 


■G 


\ m oJ 


N D =c r G c 


( 11 ) 


where c 2 = — . Values for the constants mo and m , as well as ci and C 2 are shown in Table II. 
m j 
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At the beginning of the test, the specimen is loaded initially so that the energy release rate at the 
front, Gn, reaches about 50% of Gn c corresponding to G%=0.392 kj/m 2 . The initial energy release 
rate is shown in the delamination onset plot of Figure 17 as a horizontal dashed red line. From the 
delamination onset curve, the number of cycles to delamination onset is determined as Ndso= 80, 
shown as a vertical dashed red line. The values corresponding to 40%, 30% and 20% of Gn, were 
also calculated using equation (11) and were included in the plot of Figure 17 (G 40 , Gso, G 20 - 
horizontal dashed lines; Nd 4 o, Nd 3 o, Nd 2 o - vertical dashed lines). 

5.2 Fatigue delamination growth 

The number of cycles during delamination growth can be obtained by solving equation (3) for 

N g 


N c= fdN-f^G^-da (12) 

As mentioned above, the specimen is loaded initially so that the energy release rate at the front, Gn, 
reaches about 50% of Gn c corresponding to a G/i max = Gjo=0.392 kJ/m in the current study as shown 
in the Paris Law plot of Figure 18 (solid red square). 

For practical applications, equation (3) can be replaced by an incremental equivalent expression 


Aa 
A N 


G n 


(13) 


where for the current study, increments of Aa= 0.1 mm were chosen. Starting at the initial 
delamination length, ao=25A mm, the energy release rates, G ljnax , were obtained for each increment, 
i, from the curve fit (open red circles and dashed red line) plotted in Figure 14. These energy release 
rate values were then used to obtain the increase in delamination length per cycle or growth rate 
Aa/AN from the Paris Law in Figure 18. The growth rate first increased before reaching a peak and 
then decreased with increasing a* as shown in Figure 19 ( Qso.max , dashed red line), where a* is the 
increase in delamination length. The values corresponding to 40%, 30% and 20% of Gn c were also 
determined using the same approach and were included in the plot of Figure 19 ( Qw.max , Q.w,ma X , and 
Q 2 o,max - dashed lines). For all load levels, the number of cycles during delamination growth, Ng, 
was calculated by summing the increments AN, ) 


k k -1 

N r =y AN.=y~G' n -Aa 

G / j 1 / j i,max 


i = 1 


i = 1 


(14) 


where k is the number of increments. The corresponding delamination length, a, was calculated by 
adding the incremental lengths, Aa, to the initial length, ao, 

k 

a = a 0 + a* = a 0 + ^ Aa = a 0 + k • Aa . (15) 

i - 1 


For the delamination growth phase, the increase in delamination length, a* is plotted in Figure 20 
for an increasing number of load cycles Ng. 
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The analyses were repeated for the case of applied displacements. The calculated growth rates 
Aa/AN for w 50 , max, w 40 , max, w so, max, and w 20 , max are plotted in Figure 21. After an initial increase the 
growth rates reached a flat peak and then decreased rapidly with increasing a* as shown in 
Figure 21. The number of cycles during delamination growth, Ng, and the corresponding increase in 
delamination length, a *, were also calculated using equations (14) and (15), respectively. The 
results are plotted in Figure 22. 


5.3 Combined fatigue delamination onset and growth 

For the combined case of delamination onset and growth, the total life, Nr, may be expressed as 

N t =N d +N g (16) 

where, No, is the number of cycles to delamination onset and Ng, is the number of cycles during 
delamination growth [8]. For this combined case, the increase in delamination length, a* is plotted 
in Figure 23 for an increasing number of load cycles, Nr, for all load levels ( Qso.max , Q4o,max, Q30,max, 
and Q2o,max ~ dashed lines). For the first N D cycles, the delamination length remains constant 
(horizontal dashed red line), followed by a growth section where - over Ng cycles - the delamination 
length increases following the Paris Law (dashed red line). The total life for applied cyclic 
displacements is plotted in Figure 24 for all levels (w so, max, w 40, max, woo.max, and w 20, max)- For the first 
No cycles, the delamination length remains constant (horizontal dashed red line), followed by a 
growth section where - over Ng cycles - the delamination length increases following the Paris Law 
(dashed red line). Once a delamination length is reached where the energy release rate drops below 
the assumed cutoff value, G t h (as shown in Figure 15), the delamination growth no longer follows 
the Paris Law (dashed grey line) and stops (horizontal dashed red line) as shown in Figure 24. For 
the applied cyclic displacements, w so, max, two cutoff values G,/,=0.08 kJ/m 2 and G, /,=(). 05 kJ/m 2 were 
assumed resulting in a shift of the cutoff as shown in Figure 24. 

5.4 Using the benchmark example to assess an automated analysis in a commercial FE code 

The computed energy release rate may serve as an initial verification step before the 
delamination growth analysis is initiated. The computed values at the initial crack tip in a 2D model 
or along the delamination front in a three-dimensional model should reach but not exceed the target 
value of Gnmax = 0.78 kJ/m . Once delamination growth starts in the model, the computed energy 
release rate should increase with increasing delamination length, a, as shown in Figures 14 and 15. 
The curve fits could therefore be used to check the computed energy release rate during 
delamination growth as demonstrated later. 

The growth rate Aa/AN, shown in Figures 18 and 21, may be used as a check for the correct 
implementation of the Paris Law provided this output is available. For the delamination growth 
analysis, the increase in delamination length, a*, is expected with increasing number of cycles, N, as 
shown in Figures 20 and 22. The curve fits can therefore be used as a benchmark. 

A delamination length prediction analysis that accounts for delamination fatigue onset as well 
as growth should yield results that closely resemble the plots in Figures 23 and 24. The curve fits 
can therefore be used as a benchmark as demonstrated later. 
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6. STATIC ANALYSIS BENCHMARKING 

6.1 Computed load/displacement behavior to verify correct input data and model response 

First, static analyses were performed with propagation disabled to ensure that the input was 
correct and that the models responded as expected. For this static analysis, an applied center load of 
0=2000 N was chosen for all models. To ensure that the model geometry and material input data 
for all models produced consistent results, the computed load-displacement behavior as shown in 
Figures 25 was evaluated. In comparison to the benchmark (solid grey circles and solid grey line) 
from Figure 8, which is based on analyses using enhanced plane strain elements (CPE4I), the results 
are in good agreement. The model where classical plane strain elements (CPE4) were used exhibits 
a slightly more compliant behavior (solid blue line) as shown in Figure 25. Results from models 
using enhanced and classical plane stress elements (CPS4I, CPS4) and solid elements (C3D8I) 
(solid green, light blue and black lines) are basically identical to the benchmark results. The results 
from the model with continuum shell elements (SC8R) exhibit a slightly stiffer response (solid red 
line) compared to the benchmark results. Based on the results, it was assumed that the geometry, 
maximum applied displacements and material input were defined correctly and all models 
adequately represented the benchmark case. 

6.2 Computed energy release rate distribution across the specimen width 

For all the models, the mode II strain energy release rate values were also computed for an 
applied center load 0=2000 N and plotted versus the normalized width, y/B, of the specimen as 
shown in Figure 26. The results were obtained from models shown in Figures 5 through 7. The 
results obtained from two-dimensional finite element models are plotted at the center of the 
specimen y/B =0. In comparison to the benchmark (solid grey circle) from Figure 8, which is based 
on analyses using enhanced plane strain elements (CPE4I), the results are in good agreement. Only 
the model where classical plane strain elements (CPE4) were used yielded results that were about 
1 1% higher (solid blue square) as shown in Figure 26. As expected, the mode I and mode III strain 
energy release rates were computed to be nearly zero and hence are not shown. The results from 
three-dimensional models (C3D8I, SC8R) indicate that qualitatively, the mode II strain energy 
release rate is fairly constant over almost the entire width of the specimen and only increases 
sharply near the edges. The results obtained from simulations where continuum shell elements 
(SC8R) had been used to model the specimen (see Figure 7b) were about 6.6% lower at the center 
of the specimen (solid red diamonds) compared to the results from solid elements (C3D8I) (open 
black triangles). These low results had not been observed previously when the energy release rate 
distribution across the width of a DCB had been computed [5]. Doubling the number of elements 
across the width did not alter the results (red Xs). Only when the number of elements used to model 
the thickness of the specimen was increased to three as in the original solid element model 
(Figure 6a) did the results obtained from continuum shell element models (red crosses) resemble the 
results obtained from solid models. However, if a model of the same fidelity through the thickness 
is required when using continuum shell elements, the advantage of simpler models and faster 
analysis is lost. Based on the results, it was assumed that, with respect to the computed energy 
release rates, the models adequately represented the benchmark case. 


12 



6.3 Results from automated delamination propagation analysis 

6.3.1 Computed delamination propagation for applied static displacements 

The propagation analysis was performed in two steps using the models shown in Figures 5 
through 7 for an initial delamination length, ao=25A mm. In the first step, a displacement of w=l .3 
mm was applied which nearly equaled the critical displacement, w cr it= 1.42 mm, determined earlier. 
In the second step, the applied displacement was increased to w=5.0 mm. For this second step, 
automatic incrementation was used in ABAQUS® and a small increment size (0.5% of the total 
step) was chosen at the beginning of the step. To reduce the risk of numerical instability and early 
termination of the analysis, a minimum allowed increment size (10' of the total step) was also 
chosen. The analysis was limited to 1000 increments. 

Initially, analyses were performed using two-dimensional planar models without stabilization or 
viscous regularization. Release tolerance values were varied between the default value (reltol= 0.2) 
and 0.5. Using reltol=02, 0.3 and 0.4, the load dropped at the critical point, but the displacement 
kept increasing with decreasing load as shown in Figure 27 (solid red, blue and green lines, 
respectively) where the computed resultant force (load Q) in the center of the ENF specimen is 
plotted versus the applied displacement w. Then, the analysis terminated early due to convergence 
problems. By increasing the release tolerance - as suggested in the ABAQUS® error in the message 
(.msg) file - it was possible to complete the analysis without an error message as shown in 
Figure 27, using a release tolerance value of 0.5 (solid purple line). As desired, the load dropped at 
the critical point, however, the displacement kept increasing with decreasing load. As the 
displacement continued to increase over 2.0 mm, the computed load/displacement path converged 
to the stable propagation branch of the benchmark result. For the stable path, a saw-tooth pattern is 
observed where the top results are in good agreement with the benchmark result. Based on the 
results, it was decided to introduce additional stabilization to obtain better agreement with the 
benchmark case. 

Based on problems identified during previous analyses, automated or static stabilization was not 
used in this study [5]. The results computed when contact stabilization (cs) was added are plotted in 
Figure 28. For a small stabilization factor (cs=lxl0' 6 ) and a release tolerance suggested in the 
handbook (reltol= 0.2) [13], the load decreased, and delamination propagation started shortly after 
reaching the critical point of the benchmark solution (thin solid red line). The load/displacement 
path then ran parallel to the constant displacement branch of the benchmark result and followed the 
stable propagation branch of the benchmark result. To reduce the observed overshoot, the release 
tolerance was reduced. For a stabilization factor of cs=lx 1 O' 6 and release tolerance reltol=(). 1 
(dashed blue line), the overshoot was reduced, and the computed load/displacement path then ran 
closer to the constant displacement branch of the benchmark result. Further reducing the release 
tolerance (reltol= 0.01) yielded results that were in excellent agreement with the benchmark for a 
wide range of stabilization factors (1x10' < cs <1x10' ) (solid black, light blue, orange and green 
lines - on top of each other). For the stable path, a saw-tooth pattern was observed where the top 
results were in good agreement with the benchmark result. 

When viscous regularization was added to help overcome convergence issues, a value of 0.2 
was used initially for the release tolerance as suggested in the handbook [13]. For a viscosity 
coefficient damv= 0.1 and a release tolerance of 0.2, the load and displacement overshot the critical 
point as shown in Figure 29 (solid red line). The load/displacement path then ran parallel to the 
constant displacement branch of the benchmark result and followed the stable propagation branch of 
the benchmark result. Lowering the viscosity coefficient (damv= 1x10' ) did not improve the results 
(dashed blue line). Further lowering the viscosity coefficient (damv=lxl0 4 ) reduced the overshoot 
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at the critical point, however, the computed displacement started to increase before the transition 
between the constant displacement branch and the stable propagation was reached (solid black line). 
Additional reduction in the viscosity coefficient alone (reltol= 0.2 and lxl0‘ 5 < damv <lxl0' 6 ) 
reduced the overshoot and shifted the results closer to the constant displacement branch of the 
benchmark (dashed light blue and green lines). Only simultaneously reducing the release tolerance 
{reltol= 0.1) yielded results that were in excellent agreement with the benchmark including the 
constant displacement branch and the stable propagation branch of the benchmark result (solid 
orange line). For the stable path, a saw-tooth pattern was observed where the top results were in 
good agreement with the benchmark result. Results did not change when plane stress elements 
(CPS4I) were used for the model (solid purple line). Additionally reducing the release tolerance 
{reltol= 0.01) did not yield any improvement (solid blue line) but increased the analysis time. 

An alternative way to plot the benchmark is shown in Figure 30 where the applied displacement 
w is plotted versus the increase in delamination length a*. This way of presenting the results is 
shown, since it may be of advantage for large structures where local delamination propagation may 
have little effect on the global stiffness of the structure and may therefore not be visible in a global 
load/displacement plot. However, extracting the delamination length a from the finite element 
results required more manual, time consuming post-processing of the results compared to the 
relatively simple and readily available output of nodal displacements and forces. The results plotted 
in Figure 30 are selected examples that were discussed above and were shown in the global 
load/displacement plots of Figures 27 through 29. The conclusions that can be drawn from this plot 
are identical to those discussed above. 

Results obtained from three-dimensional models are shown in Figures 31 through 33. Based on 
the results from two-dimensional planar models shown above, contact stabilization or viscous 
regularization was added for all analyses to help overcome convergence issues. Analyses where 
viscous regularization was used (lxl0' 2 < damv <1x1 0‘ 6 ) did not converge and terminated early for a 
range of release tolerance values (0.5< reltol <0.2) that had yielded converged results earlier for 
two-dimensional planar finite element models. The results computed when a small stabilization 
factor (cs=lx l0‘ 6 ) was added are plotted in Figure 31. For an initial release tolerance value 
(reltol= 0.2) as suggested in the user’s manual [13], the load dropped (solid blue line) and 
delamination propagation started shortly before reaching the critical point of the benchmark solution 
(solid grey line). The load/displacement path then ran parallel to the constant displacement branch 
of the benchmark result. At the transition between the constant displacement branch and the 
stable propagation branch of the benchmark result, the computed load was about 13% lower 
compared to the benchmark. As displacements increased, the results closely followed the stable 
propagation branch of the benchmark result. To attempt to improve the results, the release tolerance 
was varied. Neither increasing {reltol= 0.5) nor reducing (reltol= 0.1) the release tolerance had a 
significant effect on the results (red solid line and green solid line, respectively). Additionally, the 
mesh was refined, first by dividing the original element length (zla=1.0 mm) in half (A; =0.5 mm) 
and keeping the number of elements constant across the width. Second, the number of elements 
across the width was doubled resulting in the FE-model shown earlier in Figure 7a. Changing the 
mesh also did not have a significant effect on the results (black solid line and orange solid line, 
respectively). Deviation from the benchmark may be explained by the fact that the benchmark 
results were created using two-dimensional planar finite element models of the ENF specimen. As 
shown in Figure 26, the three-dimensional models yield an energy release rate distribution where 
the peak values near the edges are somewhat higher than in the center and higher than the results 
from two-dimensional planar finite element models. Therefore, delamination propagation is 
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expected to start before the peak obtained from two-dimensional planar finite element models is 
reached which shifts the entire results plot towards lower loads. 

Results obtained from models where continuum shell elements (SC8R) were used to model the 
ENF specimen (Figure 7b) are plotted in Figure 32. For an initial release tolerance value of 
reltol=0A (solid light-blue line), the load dropped and delamination propagation started 
significantly before reaching the critical point of the benchmark solution (solid grey line). In an 
attempt to improve the results, the mesh was refined by doubling the nodes in the width direction. 
Changing the mesh increased the critical load, but overall did not significantly improve the results 
(solid orange line) compared to the benchmark or to results using solid brick elements (C3D8I) as 
shown in Figure 31. Selected examples that were discussed above are plotted in Figure 33, where 
the applied displacement w is plotted versus the increase in delamination length a*. The discrepancy 
between the results obtained from models where solid brick elements (C3D8I) had been used (solid 
red, blue and green lines) and results obtained from models where continuum shell elements (SC8R) 
had been used (solid light blue and orange lines) is significant. This type of discrepancy had not 
been observed previously when propagation for a DCB specimen had been computed [3]. The 
energy release rate distributions, as shown in Figure 26, indicate that the values obtained from 
models where continuum shell elements (SC8R) had been used are lower in the center and the same 
at the edge compared to values obtained from models where solid brick elements (C3D8I) had been 
used. Based upon the values at the edge, propagation should start in continuum shell and solid 
models at the same displacement. Based upon the values in the center, propagation should start in 
the shell-element models for higher displacements. The reason for the shift of the results towards 
smaller displacements as shown in Figure 33 (solid light blue and orange lines) remains unclear and 
requires further investigation and discussion with the finite element code developers. These results 
highlight the importance of benchmarking to detect specific issues in a finite element 
implementation. 

6.3.2 Computed delamination propagation for applied static center load 

The propagation analysis was performed in two steps using the models shown in Figures 5 
through 7 for an initial delamination length, ao=25A mm. In the first step, a center load Q=\ 500 N 
was applied which equaled nearly the critical load, Q cr it= 1526 N, determined earlier. In the second 
step, the total load was increased ((7=1800 N). Automatic incrementation was used with a small 
increment size at the beginning (0.5% of the total step) and a very small minimum allowed 
increment (10‘ of the total step) to reduce the risk of numerical instability and early termination of 
the analysis. The analysis was limited to 5000 increments. 

The same steps discussed in the section on applied displacement were followed. Initially, 
analyses were performed using two-dimensional planar models without stabilization or viscous 
regularization. Release tolerance values were varied between the default value ( reltol=02 ) and 0.7. 
For the default value {reltol= 0.2), the load increase stopped at the critical point, but the analysis 
terminated immediately due to convergence problems. Even by increasing the release tolerance - as 
suggested in the ABAQUS® error in the message (.msg) file - it was not possible to complete the 
analysis without termination and reach the stable propagation path of the benchmark case. 
Therefore, additional stabilization had to be introduced in order to obtain agreement with the 
benchmark case. 

The results computed when contact stabilization (cs) was added are plotted in Figure 34. For 
small stabilization factors (cs=lx l0‘ 4 and cs=lxl0‘ 6 ) and a small release tolerance (reltol= 0.01), the 
load increased up to the critical point and delamination propagation started while the load remained 
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constant (solid black and green lines). For the stable propagation path, the results were in good 
agreement with the benchmark result. Computed results did not change when plane stress elements 
(CPS4I) were used for the model and the stabilization factor was kept constant at a low value 
(cv=lxl0" 6 ). Initially, the release tolerance value was set at the default value (reltol= 0.2) (dashed 
orange line). Then, the release tolerance was reduced to reltol= 0.1 (dashed purple line) and further 
to reltol= 0.01 (dashed green line). The results were in good agreement with the benchmark results 
over the entire load/displacement range. 

When viscous regularization was added to help overcome convergence issues, a value of 0.2 
was used initially for the release tolerance as suggested in the handbook [13] in combination with a 
viscosity coefficient damv= 0.1. The computed load/displacement behavior for this parameter 
combination was in good agreement with the benchmark case as shown in Figure 35 (solid red line). 
Reducing the viscosity coefficient to values that previously yielded good results (lxl0' 4 < damv 
<lxl0' 6 ) caused the analysis to terminate. Reducing the release tolerance to reltol= 0.1 and using 
damv= 0.1 yielded results (solid blue line) that were in excellent agreement with the benchmark 
(solid grey line). Further reducing the release tolerance to reltol= 0.01 however, caused the analysis 
to terminate. Results did not change when plane stress elements (CPS4I) were used for the model 
(reltol= 0.2 and damv=0. 1 ) as shown by the dashed red line. 

An alternative way to plot the benchmark is shown in Figure 36 where the applied center load Q 
is plotted versus the increase in delamination length a*. This way of presenting the results may be of 
advantage for large structures where local delamination propagation may have little effect on the 
global stiffness of the structure and may therefore not be visible in a global load/displacement plot. 
However, extracting the delamination length a from the finite element results required more 
manual, time consuming, post-processing of the results compared to the relatively simple and 
readily available output of nodal displacements and forces. The results plotted in Figure 36 are 
selected examples that were discussed above and were shown in the global load/displacement plots 
of Figures 34 and 35. The conclusions that can be drawn from the plot in Figure 36 are identical to 
those discussed in the two paragraphs above for Figures 34 and 35. 

Results obtained from three-dimensional models are shown in Figures 37 and 38. Based on the 
results from two-dimensional planar models shown above, only contact stabilization was added to 
help overcome convergence issues. The load/displacement result computed when a small 
stabilization factor (cs=lxl0" 6 ) was added is plotted in Figure 37. A release tolerance value 
{reltol= 0.1) was used which had yielded good results in the case where two-dimensional planar 
models had been used for the analyses. The load increase stopped (solid red line) and delamination 
propagation started shortly before reaching the critical point of the benchmark solution (solid grey 
line). As mentioned earlier, deviation from the benchmark may be explained by the fact that the 
benchmark results were created using two-dimensional planar finite element models of the ENF 
specimen. As shown in Figure 26, the three-dimensional models yield energy release rate 
distributions where the peak values near the edges are somewhat higher than in the center and 
higher than the results from two-dimensional planar finite element models. Therefore, delamination 
propagation is expected to start before the peak obtained from two-dimensional planar finite 
element models is reached which shifts the entire results plot towards lower loads. The result 
discussed was also plotted in Figure 38, where the applied center load Q is plotted versus the 
increase in delamination length a*. The conclusions that can be drawn from this plot are identical to 
those discussed above. 

In summary, good agreement between analysis results and the benchmark could be achieved for 
different release tolerance values in combination with contact stabilization or viscous 
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regularization. Selecting the appropriate input parameters, however, was not straightforward and 
often required several iterations in which the parameters had to be changed. Increasing the 
release tolerance, as suggested in the user's manual [13], may help to obtain a converged 
solution, however, leads to an undesired overshoot of the computed result compared to the 
benchmark. Therefore, a combination of release tolerance and contact stabilization or viscous 
damping is suggested to obtain more accurate results. A gradual reduction of the release 
tolerance and contact stabilization or viscous damping over several analyses is suggested. 

6.3.3 Computed delamination front shape 

Besides matching the load displacement behavior of benchmark results, a delamination 
propagation analysis should also yield a delamination front shape that is representative of the actual 
failure. An example of delamination front shapes observed by opening a tested ENF specimen are 
shown in Figure 39 [6]. Starting from the initial straight delamination front which is formed by the 
edge of the polytetrafluoroethylene (PTFE) insert, the delamination appears to grow uniformly 
across the width and remains straight during the static loading used to pre-crack the specimen [6]. 
The front observed after cyclic loading does not appear uniform across the width. The front is 
somewhat jagged, suggesting that growth happens in one location then stops and continues at 
another location across the width. On average, however, the front remains straight compared to the 
pronounced thumbnail shaped fronts observed in a Double Cantilever Beam (DCB) Specimen [3]. 

A straight front across almost the entire width of the specimen is to be expected when looking at 
the computed mode II strain energy release rate distribution plotted in Figure 26. Due to peak G- 
values near the edges, delamination propagation is expected to start in those areas first. 
Delamination propagation computed using the model with a uniform mesh across the width 
(Figure 6) is shown in Figure 40. Plotted on the bottom surface (defined in Figure 6) are the 
contours of the bond state, where the delaminated section appears in blue and the intact (bonded) 
section in red. The transition between the colors (in green/yellow) indicates the location of the 
delamination front. The initial straight front is shown in Figure 40a and was added for clarification 
in Figure 40c. The first propagation was observed at the edges of the specimen (Figure 40b) as 
expected from the distribution of the energy release rate (Figure 26). Eater, the delamination 
propagated across the width of the specimen with a jagged front shape as shown in Figure 40c. 
Propagation occurred in one location then stopped and continued at another location across the 
width. Subsequently, the entire front moved forward and, on average, remained straight. To improve 
the results, the element length, Aa, was divided in half (Act =0.5 mm). The refined model and results 
are shown in Figure 41. Note that the number of elements was kept constant across the width. As 
before, delamination propagation started at the edges of the specimen (Figure 41b). Eater, when the 
delamination propagated across the width of the specimen, the front shape still appeared somewhat 
jagged as shown in Figure 41c. To further improve the results, the number of elements across the 
width was doubled, resulting in the FE-model shown earlier in Figure 7a. As for the previous two 
models, delamination propagation started at the edges of the specimen (Figure 42a). Eater, when the 
delamination propagated across the width of the specimen, the front was jagged locally as shown in 
Figure 42b but on average moved forward as a straight front across the entire width, which is in 
agreement with the front shape observed after testing. 
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7. ASSESSMENT OF AUTOMATED DELAMINATION ONSET AND GROWTH 
ANALYSIS UNDER CYCLIC LOADING 

7.1 Static analysis to verify correct input data and model response 

First, static analyses were performed with growth disabled to ensure that the input was correct 
and that the models responded as expected. For this static analysis, only a single load cycle was 
analyzed, starting at the previously applied maximum displacement w=1.0 mm. This step was 
performed to check that the amplitude input was defined correctly and resulted in the desired 
periodic cyclic loading during the analysis. Therefore, the displacement, w, obtained as analysis 
output, was plotted versus the step time, t, as shown in Figure 43. The models using plane strain 
elements (CPE4I) (open red circles) yielded the same output as the models using plane stress 
elements (CPS4I) (blue Xs) and solid elements (C3D8I) (green crosses). Based on a comparison 
with the desired fatigue loading (grey line) it is assumed that the amplitude input was defined 
correctly and the increments are small enough to adequately represent the desired periodic fatigue 
loading shown in Figure 16. 

Second, the computed mode II energy release rate - also obtained as analysis output - was 
plotted versus the step time, t, as shown in Figure 44. As desired, the energy release rate cycled 
between the expected maximum value Gso,max = 0.5 Gu c = 0.392 kJ/m and minimum value 
Gso,min = 0.1 G so, max = 0.00392 kJ/m 2 . The model using plane strain elements (CPE4I) (open red 
circles and red line) yielded the same output as the models using plane stress elements (CPS4I) 
(blue Xs and solid blue line), which is consistent with the observations made above. The results 
from the solid model (C3D8I), taken in the center of the specimen aty=0.0 (green crosses and solid 
green line) are, as expected, in agreement with the results obtained from the two-dimensional planar 
models. 

7.2 Assessment of automated delamination onset and growth analysis 

In Figures 45 to 54, the increase in delamination length, a*, is plotted versus the total number of 
cycles, Nt, for different input parameters and models. Input parameters were varied to study the 
effect on the computed onset and growth behavior during the analysis. Based on previous results 
[5], it was decided to use an initial time increment of /o=0.001 (1/200 of a single loading cycle, 
4=0. 2s), for all analyses. Also, the solution controls were kept at fixed values. For all results shown, 
the analysis stopped when a cycle limit - used as input to terminate the analysis - was reached. 
Further details about the input parameters are discussed in the appendix where a sample input file is 
also provided. 

7.2.1 Computed delamination onset and growth for applied displacement 

Initially, analyses were performed for a benchmark case where a maximum cyclic 
displacement, W 50 , corresponding to Gso.mcvT 0.5 Gn c = 0.392 kJ/m , was applied. Two-dimensional 
planar models made of plane strain elements (CPE4I) with an element length, Aa= 1.0 mm, at the 

O 

delamination tip were used. The analysis stopped when a cycle limit, 10 , was reached. 

First, the input to define the cyclic loading was varied. In order to define the cyclic applied 
displacement, w, a set of parameters needs to be determined as shown in the equation in Figures 16 
and 43. A Fourier series is used in ABAQUS® Standard during the analysis to approximate the 
periodic cyclic loading. Based on previous results [5], it was decided to use 50 Fourier terms for 
the analysis to approximate the periodic cyclic loading. Note, that number of Fourier terms is 
required input in addition and independent of the definition of the applied cyclic displacement 
defined in equation 10 and shown in Figure 16. Further details about the required input parameters 
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are discussed in the appendix where a sample input file is also provided. The selection of the 
starting time, to, causes a phase-shift (see Figure A5 in the appendix for details). The results 
obtained for to =0.0 and to= -0.05 (open red circles and open green diamond) were in excellent 
agreement with the benchmark results (solid grey line), as shown in Figure 45. The results obtained 
for to=0.05 (open blue squares) qualitatively follows the benchmark but is shifted towards higher 
number of cycles. For all results shown, the predicted onset occurs after the benchmark result onset, 
Nd5o = 80 cycles. The threshold cutoff, where delamination growth is terminated and the 
delamination length remains constant, is predicted close to the number of cycles defined by the 
benchmark. 

Second, the effect of mesh size was studied. Therefore, the element length, Aa, at the 
delamination tip was varied as shown in the finite element models depicted in Figure 5. For this 
analysis, the cutoff value was also changed from G)/,=0.08 kJ/m to G,/,=0.05 kJ/m (see Figure 24). 
This change served as an additional checkpoint to assure that the automated growth was terminated 
and the delamination length remained constant once the cutoff value was reached. The results 
obtained from the respective models are plotted in Figures 46. Excellent agreement with the 
benchmark curve (grey solid line) could be achieved for all element lengths studied. For smaller 
element length, Aa=025 mm (open green diamonds and solid green line), the predicted onset occurs 
at a slightly lower number of cycles compared to the larger element length of Aa= 1.0 mm (open red 
circles and solid red line) and Aa= 0.5 mm (open blue squares and solid blue line). The observed 
mismatch is largely due to the increased element length, which causes the first growth step to be 
larger compared to the results obtained from shorter elements. For all models, the threshold cutoff, 
where delamination growth is expected to stop and the delamination length remains constant, is 
predicted close to the number of cycles defined by the benchmark. 

The analyses were repeated for the other benchmark cases with applied maximum cyclic 
displacement, W 40 (dashed red line), W 30 (dashed blue line), and W 20 (dashed green line), as shown in 
Figure 47. Two-dimensional planar models - made of plane strain elements (CPE4I) with an 
element length, Aa= 1.0 mm, at the delamination tip - were used. As before, it was decided to use 50 
Fourier terms to approximate the periodic cyclic loading. A starting time, to= 0.0, was selected 
based on the results from above. The computed results for W 40 (open red circles and solid red line), 
W30 (open blue squares and solid blue line), and W20 (open green diamond and solid green line) are in 
excellent agreement with the respective benchmark cases as shown in Figure 47. 

The results obtained for models made of three-dimensional solid elements (shown in Figures 6 
and 7) are presented in Figure 48. The observed delamination fronts during growth, which will be 
discussed later, were not as uniform across the width of the specimen as would be expected from the 
static propagation (see Figures 38-42) or from the distribution of the energy release rate across the 
width (see Figure 26). Therefore, the increase in delamination length, a* was plotted for three 
selected locations: the center of the specimen at y=0.0 (open circles); edge 1 at y=0.5 (open 
squares); and edge 2 at y= -0.5 (open diamonds), as shown in Figure 48. To reduce the 
computational effort, the analysis was initially performed for an element length, da=1.0 mm and 
only one Fourier term to approximate the periodic cyclic loading (results in red - CPU time 8.23 
10 5 s or 9.5 days ). Good agreement with the benchmark curve (grey solid line), however, could be 
achieved only when the number of Fourier terms to approximate the periodic cyclic loading was 


The analyses were performed on a Quad-Core AMD Opteron(tm) Processor 8378 running openSUSE 1 1.3 (x86_64) 
using ABAQUS® Standard 6.9 and 6. 10 on a single CPU. 
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increased to 50 as before (results in blue), which increased the computation time somewhat to 8.60 
10 5 s (4.3 % increase). 

To improve the results, the element length, Aa, was divided in half (Aa= 0.5 mm). Note that the 
number of elements was kept constant across the width. The results obtained when one Fourier term 
was used deviated from the benchmark as shown in Figure 49 (in green). Increasing the number of 
Fourier terms to approximate the periodic cyclic loading to 50 improved the results somewhat 
(purple lines and open symbols). The onset and initial growth were captured sufficiently accurately. 
Results for all locations (center and both edges), however, started to deviate from the benchmark 
after the initial growth phase of the analysis. To further improve the results, the number of elements 
across the width was doubled resulting in the FE-model shown earlier in Figure 7a. The results 
started to deviate from the benchmark very early in the analysis (in black). The observed 
performance for the refined meshes remains unclear and further discussion with the software 
developers is required. 

7.2.2 Computed delamination onset and growth for applied center load 

Analyses were also performed for benchmark cases where a maximum cyclic center load, Q, 
was applied to the model of the specimen. Two-dimensional planar models made of plane strain 
elements (CPE4I) with an element length, Aa= 1.0 mm, at the delamination tip were used first and a 
center load, Q 50 , corresponding to Gso.mcoT 0.5 Gn c = 0.392 kJ/m was applied. The analysis stopped 
when a cycle limit of 6000 was reached. 

First, the input to define the cyclic loading was varied. In order to define the cyclic applied 
load, Q, a set of parameters needs to be determined as shown in the equation in Figures 16 and 43. 
As mentioned above, it was decided to use 50 Fourier terms to approximate the periodic cyclic 
loading. The selection of the starting time, to, causes a phase-shift (see Figure A5 in the appendix 
for details). The results obtained for to—0.05 and to=0.0 (open red circles and open green diamond) 
were in excellent agreement with the benchmark results (solid grey line), as shown in Figure 50. 
The results obtained for ^=0.05 (open blue squares) qualitatively follows the benchmark but is 
shifted towards a higher number of cycles. For all results shown, the predicted onset occurs after the 
benchmark result onset, Adjo= 80 cycles. 

Second, the effect of mesh size was studied. Therefore, the element length, Aa, at the 
delamination tip was varied as shown in the finite element models depicted in Figure 5. The results 
obtained from the respective models are plotted in Figure 51. Excellent agreement with the 
benchmark curve (grey solid line) could be achieved for both element lengths studied. For the 
smaller element length, Aa= 0.5 mm, (open blue squares and solid blue line), the predicted onset 
occurs at a slightly lower number of cycles compared to an element length of Aa= 1.0 mm (open red 
circles and solid red line) and is closer to the benchmark. The observed mismatch is largely due to 
the increased element length, which causes the first growth step to be larger compared to the results 
obtained from shorter elements. 

The analyses were repeated for the other benchmark cases with applied maximum cyclic center 
loads, Q 40 (dashed red line), Q 30 (dashed blue line), and Q 20 (dashed green line), as shown in 
Figure 52. Two-dimensional planar models - made of plane strain elements (CPE4I) with an 
element length, Aa= 1.0 mm, at the delamination tip - were used. As before, it was decided to use 50 
Fourier terms to approximate the periodic cyclic loading. A starting time, to= 0.0, was selected 
based on the results from above. The computed results for Q 40 (open red circles and solid red line), 
Qso (open blue squares and solid blue line), and Q 20 (open green diamond and solid green line) are 
in excellent agreement with the respective benchmark cases as shown in Figure 52. 
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The results obtained for models made of three-dimensional solid elements (shown in Figures 6 
and 7) are plotted in Figure 53. The observed delamination fronts during growth, which will be 
discussed later, were not as uniform across the width of the specimen as would be expected from the 
static propagation (see Figures 38-42) and the distribution of the energy release rate across the width 
(see Figure 26). Therefore, the increase in delamination length, a* was plotted for three locations: 
the center of the specimen at y=0.0 (open circles); edge 1 at y= 0.5 (open squares); and edge 2 at 
y= -0.5 (open diamonds), as shown in Figure 53. As before, the analysis was initially performed for 
an element length, da=1.0 mm and only one Fourier term to approximate the periodic cyclic 
loading in order to reduce the computational effort (results in red - CPU time 2.8 10 5 s). The results 
in the center and at edge 1 (y= 0.5) were in good agreement with the benchmark curve (grey solid 
line). Upon increasing the number of Fourier terms to 50, the results in the center were not 
significantly affected (open circles). However, the results at the edges reversed with edge 2 (y= -0.5) 
edge being in good agreement with the benchmark curve (open blue diamonds) and the results at 
edge 1 (y= 0.5) slowing prematurely (open blue squares). The observed performance remains 
unclear and further discussion with the software developers is required. The corresponding 
delamination front shapes during growth will be discussed later. 

To improve the results, the element length, Aa, was divided in half (Aa= 0.5 mm). Note that the 
number of elements was kept constant across the width. The results obtained when one Fourier term 
was used deviated from the benchmark as shown in Figure 54 (in green). The results in the center 
and at the edges were shifted towards a lower number of cycles compared to previous results from 
two-dimensional analyses and the benchmark (solid grey line). Increasing the number of Fourier 
terms to 50 did not improve the results (purple line and open symbols). Only the onset and initial 
growth was captured sufficiently accurately. The computational effort was substantial and reached 
1.83 10 6 s (21 days). For further improvement, the number of elements across the width was 
doubled resulting in the FE-model shown earlier in Figure 7a. This model had yielded results that 
had been in excellent agreement with the benchmark in the static case discussed above. However, 
for the analysis of onset and growth under a cyclic center load the results only followed the 
benchmark during the onset and initial growth phase (black lines and open symbols) and then 
quickly deviated from the benchmark. The analysis was terminated by the user after 1.44 1 0 6 s (16.6 
days). The corresponding delamination fronts during growth will be discussed later. The observed 
performance remains unclear and further discussion with the software developers is required. Long 
analysis times may be avoided if the node release is made independent of the complicated cyclic 
fatigue implementation in ABAQUS® Standard. 

7.2.3 Computed delamination front shape 

Besides matching the delamination length versus number of cycles behavior of the benchmark 
results, a delamination onset and growth analysis should also yield a delamination front shape that is 
representative of the actual failure observed during testing. An example of delamination front 
shapes observed by opening a tested ENF specimen was shown previously in Figure 39 [6]. The 
front observed after cyclic loading does not appear uniform across the width. The front is somewhat 
jagged, suggesting that growth happens in one location, then stops and continues at another location 
across the width. On average, however, the front remains straight compared to the pronounced 
thumbnail shaped fronts observed in a Double Cantilever Beam (DCB) Specimen [3]. A straight 
front across almost the entire width of the specimen is to be expected when looking at the computed 
mode II strain energy release rate distribution plotted in Figure 26. Due to peaks near the edges 
delamination propagation is expected to start in those areas. 
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Delamination growth computed using the model with a uniform mesh across the width 
(Figure 6), applied displacement, W 50 , and one Fourier term used in the analysis is shown in 
Figure 55. Plotted on the bottom surface (defined in Figure 6) are the contours of the bond state, 
where the delaminated section appears in blue and the intact (bonded) section in red. The transition 
between the colors (in green/yellow) indicates the location of the delamination front. The initial 
straight front was added for clarification in Figures 55a, c, d and e. The first growth was observed at 
the edges of the specimen (Figure 55b) as expected from the distribution of the energy release rate 
(Figure 26). The front shown in Figure 55b corresponds to a total number of Nr =400 cycles. Later, 
the delamination propagated across the entire width of the specimen and the front started to advance 
in the center first as shown in Figure 55c. Further growth occurred at the center and edge 1 aty=0.5, 
which resulted in a curved front as shown in Figure 55d. Later, growth occurred predominately at 
edge 2 (y=- 0.5), resulting in the slightly slanted final front after A^r = 1 0 8 cycles shown in 
Figure 55e. The images of the delamination fronts help visualize the results plotted previously in 
Figure 48. 

To improve the results, 50 Fourier terms were used as discussed earlier and the resulting 
computed growth is shown in Figure 56. As before, the first growth was observed at the edges of the 
specimen (Figure 56b) as expected from the distribution of the energy release rate (Figure 26). 
Later, the delamination propagated across the entire width of the specimen and the front started to 
advance in the center first as shown in Figure 56c. As the delamination grew, the front always 
turned into a straight front across the entire width before new growth started from the center as 
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shown in Figure 56d. The final front after Nt =10 cycles shown in Figure 56e is in good agreement 
with the experimental finding shown in Figure 39. 

For a smaller element length ( Act =0.5 mm), selected details of the growing fronts are shown in 
Figures 57 through 59. Delamination growth computed using applied displacement, W 50 , and one 
Fourier term used in the analysis is shown in Figure 57. Note that the number of elements was kept 
constant across the width. As before, delamination growth initially started at the edges of the 
specimen. Later, the delamination grew across the width of the specimen, however, the center and 
both edges lagged behind as shown in Figure 57a. The front retained this double-curved shape 
throughout the analysis, as shown in Figure 57b. To improve the results, 50 Fourier terms were used 
in the analysis and the resulting fronts are shown in Figure 58. Again, delamination growth initially 
started at the edges of the specimen. As the delamination grew across the width of the specimen, 
however, the center and edge 1 (y=0.5) lagged behind as shown in Figure 58a. Later, growth 
occurred predominately at that edge, resulting in the changed curved front shown in Figure 58b. To 
further improve the results, the number of elements across the width was doubled, resulting in the 
FE-model shown earlier in Figure 7a. As for the previous models, delamination propagation started 
initially at the edges of the specimen. Later, when the delamination propagated across the width of 
the specimen the front is jagged locally as shown in Figure 59a. The front remained irregular as 
shown in Figure 59b until the analysis was terminated after 35 days. The observed performance for 
the refined meshes remains unclear and further discussion with the software developers is required. 

Delamination growth computed using the model with a uniform mesh across the width (A/ =1.0 
mm, Figure 6); applied load, Q 50 ', and one Fourier term used in the analysis is shown in Figure 60. 
The initial straight front was added for clarification in Figures 60a, c, d and e. The first growth was 
observed at the edges of the specimen (Figure 60b) as expected from the distribution of the energy 
release rate (Figure 26). Later, the delamination propagated across the entire width of the specimen 
and the front started to advance faster at edge 2 (y= -0.5) developing into a somewhat slanted front, 
as shown in Figure 60c. Then growth occurred more in the center and at edge 1 (y=0.5), which 
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resulted in a curved front, as shown in Figure 60d. The final front after N T = 6000 cycles shown in 
Figure 60e remained curved with the growth at edge 2 (y= -0.5) lagging behind. The images of the 
delamination fronts help visualize the results plotted in Figure 53. 

For a smaller element length ( Act =0.5 mm), selected details of the growing fronts are shown in 
Figures 61 through 63. Delamination growth computed using applied load, Qw, and one Fourier 
term used in the analysis is shown in Figure 61. Note that the number of elements was kept constant 
across the width. As before, delamination growth initially started at the edges of the specimen. 
Later, the delamination started to advance faster at edge 1 (y=0.5), developing into a somewhat 
slanted front as shown in Figure 61a. The front shape later advanced at edge 2 (y= -0.5) and the 
center and turned into a double-curved front as shown in Figure 61b at the end of the analysis after 
6000 cycles. To improve the results, 50 Fourier terms were used in the analysis and the resulting 
fronts are shown in Figure 62. Again, delamination growth initially started at the edges of the 
specimen. As the delamination grew across the width of the specimen, however, the delamination 
started to advance faster at edge 1 (y=0.5), developing into a somewhat slanted front, as shown in 
Figure 62a. Later, growth occurred predominately at edge 2 (y= -0.5), resulting in the final curved 
front after Nr =6000 cycles shown in Figure 62b. To improve the results, the number of elements 
across the width was doubled resulting in the FE-model shown earlier in Figure 7a. As for the 
previous models, delamination propagation started initially at the edges of the specimen as shown in 
Figure 63a. Later, when the delamination propagated across the width of the specimen, the front 
was jagged locally and remained irregular as shown in Figure 63b until the analysis was terminated 
after 1.44 1 0 6 s CPU time (Nr= 9540 cycles). The observed performance for the refined meshes 
remains unclear and further discussion with the software developers is required. 


8. SUMMARY AND CONCLUSIONS 

The development of benchmark examples, which allow the assessment of the static 
delamination propagation as well as the growth prediction capabilities was presented and 
demonstrated for ABAQUS® Standard. The example is based on a finite element model of mode II 
End Notched Flexure (ENF) specimen. The model is independent of the analysis software used and 
allows the assessment of the delamination growth prediction capabilities in commercial finite 
element codes. 

First, static benchmark results were created based on the approach developed in reference 3, 
using two-dimensional finite element models for simulating the ENF specimens with different 
initial delamination lengths, ao. For each delamination length modeled, the load and displacements 
were monitored. The mode II strain energy release rate was calculated for a fixed applied load. It 
was assumed that the delamination propagated when the mode II strain energy release rate reached 
the fracture toughness value. Thus, critical loads and critical displacements for delamination 
propagation were calculated for each initial delamination length modeled. From these critical 
load/displacement results, benchmark solutions were created. It was assumed that the 
load/displacement relationship computed during automatic propagation should closely match the 
benchmark cases. 

Second, the development of a benchmark example for delamination fatigue growth prediction 
was presented step by step. The number of cycles to delamination onset was calculated from the 
material data for mode II fatigue delamination growth onset from reference 6. The number of cycles 
during stable delamination growth was obtained incrementally from the material data for mode II 
fatigue delamination propagation. For the combined benchmark case of delamination onset and 
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growth, the delamination length was calculated for an increasing total number of load cycles. 
Second, starting from an initially straight front, the delamination was allowed to grow under cyclic 
loading. The number of cycles to delamination onset and the number of cycles during stable 
delamination growth for each growth increment were obtained from the analysis. 

After creating the benchmark cases, the approach was demonstrated for the commercial finite 
element code ABAQUS®. Starting from an initially straight front, the delamination was allowed to 
propagate under static loading or grow under cyclic loading based on the algorithms implemented 
into the software. Input control parameters were varied to study the effect on the computed 
delamination propagation and growth. 

The results showed the following: 

• The benchmarking procedure was capable of highlighting the issues associated with 
the input parameters of a particular implementation. 

• In general, good agreement between the results obtained from the propagation and 
growth analysis and the benchmark results could be achieved by selecting the 
appropriate input parameters. However, selecting the appropriate input parameters 
was not straightforward and often required an iterative procedure. 

• Different sets of input parameters were identified to be important for the automated 
propagation analysis under static loading compared to the automated onset and 
growth analysis under cyclic loading. 

• The results for automated delamination propagation analysis under static loading 
showed the following: 

o Increasing the release tolerance as suggested in the ABAQUS® handbook 
helped to obtain a converged solution, however, it led to undesired overshoot 
of the result. 

o A combination of release tolerance and contact stabilization or viscous 
damping is required to obtain more accurate results. 

o A gradual reduction of the release tolerance and contact stabilization or 
viscous damping over several analyses is suggested. 

o Results obtained from continuum shell models were offset compared to results 
obtained from three-dimensional solid models. Energy release rates were lower 
in the center of the specimen. Values were higher at the edges causing 
propagation to begin at lower loads and displacements than the solid models or 
the benchmark case. 

o Computed delamination fronts obtained from three-dimensional solid models 
were somewhat jagged. Fronts obtained from refined meshes matched the 
experimentally observed shapes. 

• The results for automated delamination growth analysis under static loading showed 
the following: 

o Stabilization or viscous damping is not required. The release tolerance has no 
effect. 

o The onset prediction appeared much more sensitive to the input parameters 
than the growth prediction. 

o Best results were obtained when a sine curve representation of the cyclic applied 
displacement was selected in combination with the starting time, to= 0.0. 
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o Consistent results were obtained when input parameters were selected such 
that 50 terms in the Fourier series were used during the execution of 
ABAQUS® Standard to approximate the periodic cyclic loading, 
o Results obtained from three-dimensional solid models appeared mesh 
dependent. Results obtained from refined meshes deviated from the benchmark 
compared to results from coarser meshes. Computed delamination fronts did not 
match the experimentally observed shape. Since the reason for this discrepancy 
is unclear, further discussion with the software developers is required, 
o The analyses for three-dimensional models of a simple ENF specimen required 
many days of computation time. Improvements to the implementation by 
making the node release independent of the complicated cyclic fatigue 
implementation may alleviate this problem. 

Overall, the results are promising and the current findings concur with previously published 
conclusions [3,4,5]. Further studies, however, are required which should include the assessment of 
the propagation capabilities in more complex mixed-mode specimens and on a structural level. 

Assessing the implementation in one particular finite element code illustrated the value of 
establishing benchmark solutions since each code requires specific input parameters unique to its 
implementation. Once the parameters have been identified, they may then be used as a starting point 
to model delamination growth for more complex configurations. 
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TABLE I. MATERIAL PROPERTIES [8]. 


Unidirectional Graphite/Epoxy Prepreg 

En = 161 GPa 

£22= H-38 GPa 

£33 = 11.38 GPa 

vi 2 = 0.32 

V13 = 0.32 

V23 = 0.45 

Gi2 = 5.2 GPa 

G13 = 5.2 GPa 

G23 = 3.9 GPa 


The material properties are given with reference to the ply coordinate axes where index 1 1 denotes the ply principal 
axis that coincides with the direction of maximum in-plane Young’s modulus (fiber direction). Index 22 denotes the 
direction transverse to the fiber in the plane of the lamina and index 33 the direction perpendicular to the plane of the 
lamina. 


TABLE II. FRACTURE PARAMETERS. 

Fracture Toughness Data [6,9] - Figures 2, A1 

Gj c = 0.21 kJ/m 2 Gjjc = G/// c =0.78 kJ/m 2 77= 2.57 

Delamination Growth Onset Data [6] - Figures 3, 17, A3, 

mo=0.78 m/=-0.16 

£7=0.213 c 2 =-6.329 

Delamination Growth Rate Data (Paris Law) [6] - Figures 4, 18, A5 

Gllc = 0.78 kJ/m 2 G t h = 0.08 kJ/m 2 

c=0.33 77=5.55 
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APPENDIX 

Delamination fatigue growth analysis in ABAQUS® 

Delamination growth at the interfaces in laminated composites subjected to cyclic loadings 
can be simulated in ABAQUS® by specifying the low-cycle fatigue criterion propagation analysis 
using the direct cyclic approach [13]. The interface along which the delamination (or crack) 
propagates must be indicated in the model using a fracture criterion definition. The onset and 
growth of fatigue delamination at the interfaces are characterized by the relative fracture energy 
release rate. The fracture energy release rates at the crack tips in the interface elements are 
calculated based on the virtual crack closure technique (YCCT). The low-cycle fatigue analysis 
in ABAQUS® is a quasi-static analysis on a structure subjected to sub-critical cyclic loading. The 
low-cycle fatigue analysis in ABAQUS® uses the direct cyclic approach to obtain the stabilized 
cyclic response of the structure directly. The direct cyclic analysis uses a combination of Fourier 
series and time integration of the nonlinear material behavior to obtain the stabilized cyclic 
response of the structure iteratively and therefore avoids the numerical expense associated with a 
transient analysis. The direct cyclic analysis in ABAQUS® is therefore suited for very large 
problems in which many load cycles must be applied to obtain the stabilized response. The direct 
cyclic analysis in ABAQUS®, however, is limited to geometrically linear behavior and fixed 
contact conditions. The theory and algorithm to obtain a stabilized response using the direct 
cyclic approach are described in detail in the ABAQUS® Theory Manual [14]. 

Required input for ABAQUS® 

The input required to perform a delamination onset and growth analysis in ABAQUS® 
Standard is discussed in the following paragraphs. It is assumed that the reader is familiar with 
ABAQUS® Standard and the syntax used in the input file (.inp). The focus is therefore on the 
specific input that relates to delamination propagation, low-cycle fatigue and the direct cyclic 
approach in ABAQUS®. An example input file is given at the end of this appendix to provide an 
overview of an entire analysis and assist the readers in creating their own analyses. 

Input for delamination propagation 

The interface along which the delamination (or crack) propagates must be indicated in the 
model using a fracture criterion definition: 

* DEBOND , SLAVE=VCCT_TOP ,MASTER=VCCT_BOT , FREQ=1 

* FRACTURE CRITERION, TYPE=f atigue , MIXED MODE BEHAVIOR=BK, TOLERANCE=<tol> 

<cl>,<c2>,<c3>,<c4>,<rl>,<r2>,<GIc>,<GIIc>, 

<GIIIc>,<eta> 

where vcct_top and vcct_bot are interface surfaces as shown in Figure A1 and <toi> is the 
tolerance within which the crack propagation criterion must be satisfied. The input parameters 
for the fracture criterion are obtained from the static mixed-mode failure criterion, the 
delamination growth onset criterion and the growth rate shown in Figures A2 -A5. 

The critical energy release rates <gic>,<giic>,<giiic> and the curve fit parameter <eta> are 
obtained from the mixed-mode failure criterion as shown in Figure A2. A quasi static mixed-mode 
fracture criterion is determined for a material by plotting the interlaminar fracture toughness, G c , 
versus the mixed-mode ratio, Gu/Gt as discussed earlier and shown in Figure A2. 

The parameters <ci>,<c2> are obtained by solving the law for growth onset (shown in Figure 
A3) for the number of cycles N as shown in equation (8). The parameters <c3>,<c4> are obtained 
directly from the Paris Law as shown in Figure A4. The parameter <ri> is calculated from the 
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energy release rate cutoff value, G t i h and the fracture toughness, Gi c , as shown in Figure A4. To 
calculate the parameter <r2>, the user needs to define an energy release rate upper limit, G p j, 
above which the fatigue crack will grow at an accelerated rate as shown in Figure A4. For the 
current benchmark example, G p i was chosen to be 90% of the fracture toughness. 

Input for cyclic loading 

Defining a low-cycle fatigue analysis using the direct cyclic approach in ABAQUS® 
Standard requires the definition of an amplitude curve which describes the relative load 
magnitude: 

* amplitude , name=test , DEFINITION=PERIODIC 

<n> , <omega> , <t 0 > , <A 0 > 

<A 1 >,<B 1 > 

where test is the label to be used to refer to the amplitude curve. The parameters defined in the 
first line are the number of terms in the Fourier series, <n>, the circular frequency, <omega>, in 
radians per time, the starting time, <t 0 >, and the constant term in the Fourier series, <Ao>, as 
shown in Figure A5. Using different starting times, <t 0 >, causes a phase-shift in the amplitude 
curve (dashed lines) as shown in Figure A5. The parameters defined in the second line are the 
first coefficient of the cosine terms, <a x >, and the first coefficient of the sine terms, <b x >, also 
shown in Figure A5. 

The amplitude curve is then referenced in the definition of the cyclic loading. In the current 
example, prescribed displacements were used to simulate the displacements of the ENF specimen: 

♦BOUNDARY, AMPLITUDE=test 

LFRONT , 1, 1, -1.0 

where lfront is the node set located at the center of the specimen where the displacements are 
applied as shown in Figure A1 . The factor l . o is used to multiply the relative magnitude defined by 
the amplitude curve (shown in Figure A5) and obtain the applied cyclic displacement, w, as shown 
in Figure A6. 

The direct cyclic approach in ABAQUS® Standard is used to obtain the stabilized cyclic 
response of a structure directly: 

♦direct cyclic, fatigue 

, , , <n i > , <n max > , <An> , <i max >, 

, , <N t > , , 

where the parameter fatigue is used to perform a low-cycle fatigue analysis. The parameters 
defined in the first line are the initial time increment, <i 0 >, the time of a single loading cycle, 
<t s > (as shown in Figure A6), the minimum time increment allowed (not used), the maximum 
time increment allowed (not used), the initial number of terms in the Fourier series, <ni>, the 
maximum number of terms in the Fourier series, <n max >, the increment in number of terms in the 
Fourier series, <An> , and the maximum number of iterations allowed in a step, <i max >. The 
parameters defined in the second line are the minimum increment in number of cycles over 
which the damage is extrapolated forward (default used), the maximum increment in number of 
cycles over which the damage is extrapolated forward (default used), the total number of cycles 
allowed in a step, <n t >, and the damage extrapolation tolerance (default used). The time of a 
single loading cycle was kept constant at 4=0.2 s for all analyses. Most analyses were run up to 

O 

44=1 0 cycles in order to reach the threshold after which delamination growth stops as shown in 
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Figure 24. All other input parameters were varied and the effect on the results was studied as 
discussed in the main part of this report. 


Control parameters direct cyclic analysis 

Solution controls in ABAQUS® Standard can be reset and modified by using keyword 

*controls,type=direct cyclic 

<I PI >, <CR n > , <CU n > , <CR 0 > , <CU 0 > , 

where the parameter direct cyclic is used to set parameters that will be used to control the 
stabilized state and plastic ratcheting detections and to specify when to impose the periodicity 
condition for direct cyclic analysis. If this keyword is omitted, default parameters are used. The 
parameters defined in the first line are the iteration number at which the periodicity condition is 
first imposed, <i PI >, (default used), the stabilized state detection criterion for the ratio of the 
largest residual coefficient on any terms in the Fourier series to the corresponding average flux 
norm, <CR n >, the stabilized state detection criterion for the ratio of the largest correction to the 
displacement coefficient on any terms in the Fourier series to the largest displacement 
coefficient, <cu n >, plastic ratcheting detection criterion for the ratio of the largest residual 
coefficient on the constant term in the Fourier series to the corresponding average flux norm, 
<CRo>, and the plastic ratcheting detection criterion for the ratio of the largest correction to the 
displacement coefficient on the constant term in the Fourier series to the largest displacement 
coefficient, <cu 0 > [13]. 


Example input files 

An input file is given to provide an overview of an entire analysis and assist the readers in 
creating their own analyses. The analysis was divided into two steps. In the first step, a small 
static preload step was introduced as a work around to avoid problems discovered with the initial 
contact conditions (ABAQUS® bug v68_1987). The second step performed the desired cyclic 
analysis. It was found that the prescribed displacements in the static preload step had to be small 
(0.001 mm) compared to the prescribed displacements ( v min =0.1, v max =-\.0), which were used to 
simulate the cyclic opening of the arms of the ENF specimen. In a previous study, analyses using 
larger preload steps did not converge in the second step and terminated prematurely [5]. 

For all analyses, the input to define the fracture criterion (<Gic>,<Giic>,<Giiic>,<eta>), the 
parameters for delamination onset (<ci>,<c2>), and delamination growth (Paris Law) 
(<c3>,<c4>,<ri>,<r2>) were kept constant. The parameters to define the load frequency 
(<omega> , <a 0 >) as well as the minimum and maximum applied displacement (<t s >, v max =- 1.0) 
were also kept constant during all analyses. Based on previous results [5], it was decided to use an 
initial time increment of io=0.00 1 (1/200 of a single loading cycle, / S =0.2s), for all analyses. Also 
the solution controls were kept at fixed values. Other parameters required to define the Fourier 
series, which is used to define the cyclic load, were varied. The ABAQUS® keywords shown in 
bold type were discussed in detail in the previous paragraphs. 


29 



Input file for fatigue onset and growth analysis 

*HEADING 

ENF-UD-IM7/8552, a=10 mm 
units: mm, N, MPa 

*** elements , nodes, material , etc 
*ELEMENT, TYPE=CPE4I , ELSET=EALL 


*NSET, NSET=BONDED, GENERATE 
205, 845, 8 

**** surface and contact definition **** 

* SURFACE, TYPE=ELEMENT, NAME=VCCT_BOT 
EL_BOT, S3 

* SURFACE, TYPE=ELEMENT, NAME=VCCT_TOP 
EL_TOP, SI 

* INITIAL CONDITION, TYPE=CONTACT 
VCCT_TOP, VCCT_BOT, BONDED 

* CONTACT PAIR, INTERACTION=VCCT, AD JUST=BONDED 
VCCT_TOP, VCCT_BOT 

* * * * VCCT STAGE I **** 

* SURFACE INTERACTION, NAME = VCCT 
<width> 

**** vcCT fatigue input 
^parameter 

** Damage and tolerance parameters 
tol=0 . 001 

** Fracture toughness: 

GIc = 0.208 
GIIc = 0.78 
GIIIc = 0.78 
** B-K parameter: 
eta=2 . 5713 

** width in the plane stress/strain direction 
width =25.4 

** fatigue crack growth data ** 
cl=0 . 213 
c2=-6 . 329 
c3=0. 33185 
c4=5 . 5519 
rl=0 . 102 
r2=0 . 9 
*** 

*** amplitude **** 

* amplitude , name=tes t , DEFINITION=PERIODIC 

1,31.42,0.05,0.55 

0,0.45 

*** history data *** 

*** static ramp up step **** 

*STEP, NLGEOM, INC= 10000 

* STATIC 
0.001, 0.001 

* * 

* DEBOND , SLAVE=VCCT_TOP ,MASTER=VCCT_BOT , FREQ=1 
* FRACTURE CRITERION , TYPE=VCCT , MIXED MODE BEHAVI OR=BK 
1 . 0e6 , 1 . 0e6 , 1 . 0e6 ,<eta> 
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* * 


**** ENF loading 

*** center deflection for ENF, 2D- full model 
* BOUNDARY , TYPE=D IS PLACEMENT 
LFRONT , 1, 1, -0.001 

* * 

*** center load for ENF, 2D- full model 
****CLOAD 

*** 408,1 , -1.0 
* * 

** field and history output ** 

*0UTPUT, FIELD, VARIABLE=PRESELECT , FREQ=1 
*Output , history, VARIABLE=PRESELECT, freq=l 
*N0DE output, NSET=LFRONT 
RF1 

*N0DE output, NSET=LFRONT 
U1 

*N0DE output, NSET=SUPPORT 
RF1 

*END STEP 

**** cyclic loading step **** 

*STEP , INC= 10000 
* direct cyclic , fatigue 
0.001,0.2,, , 50, 50, , 20 
3,6,100000000, , 

* DEBOND , SLAVE=VCCT_TOP ,MASTER=VCCT_BOT , FREQ=1 

* FRACTURE CRITERION , TYPE=fatigue , MIXED MODE BEHAVIOR=BK , TOLERANCE=<tol> 
<cl> , <c2> , <c3> , <c4>, <rl> , <r2> , <GIc> , <GI Ic> , 

<GIIIc>,<eta> 

*** run analysis first with default values 

^controls , type=direct cyclic 

, 100 , IE-3 , 100 , IE-3 

**** ENF loading 

^BOUNDARY , AMPLITUDE=test 

LFRONT, 1, 1, -1.0 

* * 

*** center load for ENF, 2D-full model 

* * *CLOAD, AMPLITUDE=tes t 

** 408,1, -682.5 

* * 

** field and history output ** 

*OUTPUT , FIELD, VARIABLE=PRESELECT , FREQ=25 
*ELEMENT OUTPUT 
cycleini, status, sdeg 

* CONTACT OUTPUT, MASTER=VCCT_BOT , SLAVE=VCCT_TOP 
dbt, dbsf , dbs , openbc, crsts , enrrt, ef enrrtr, bdstat 
*Output, history, VARIABLE=PRESELECT, freq=25 
*NODE output, NSET=LFRONT 
RF1 

*NODE output, NSET=LFRONT 
U1 

*NODE output, NSET=SUPPORT 
RF1 

*END STEP 
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dimensions 



Figure 1. End-Notched Flexure Specimen (ENF) [6], 


G , 
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kJ/m 2 



mixed mode ratio G /G 

Figure 2. Mixed mode fracture criterion. 
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Figure 3. Delamination growth onset for mode II for IM7/8552 [6], 
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Figure 4. Delamination growth rate (Paris Law) for mode II for IM7/8552 [6], 



o experimental data [6] 
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a. Deformed FE-model ofENF specimen with initial delamination before growth. 




Aa=0.5 mm 


c. Delamination front detail of a refined FE-model of an ENF specimen. 



Aa=0.25 mm 


d. Delamination front detail of a further refined FE-model of an ENF specimen. 
Figure 5. Two-dimensional finite element model of an ENF specimen. 
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a. Deformed three-dimensional FE-model with initial delamination front before growth. 



b. Detail of three-dimensional FE-model around delamination front. 
Figure 6. Full three-dimensional finite element model of an ENF specimen. 
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a. Deformed three-dimensional model of an ENF specimen with a fine mesh before propagation. 



b. Deformed continuum shell model of an ENF specimen after delamination propagation. 
Figure 7. Three-dimensional finite element models of an ENF specimen. 
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Figure 8. Load-displacement behavior for ENF specimens with 
different initial delamination lengths, a . 
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Figure 9. Calculated critical load-displacement behavior for ENF specimen. 
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Figure 10. Benchmark cases for applied load, Q, and displacement, w,forENF specimen. 



Figure 11. Benchmark case for applied displacement, w, plotted versus 
increase in delamination length, a*. 
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Figure 12. Benchmark case for applied load, Q, plotted versus increase in delamination length, a * 



displacement w, mm 


Figure 13. Maxium cyclic loads and applied displacments for an ENF specimen. 
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Figure 14. Energy release rate - delamination length behavior 
for ENF specimen with constant load. 
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Figure 15. Energy release rate - delamination length behavior 
for ENF specimen with constant displacement. 
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Figure 16. Applied cyclic center deflection for ENF specimen. 



Figure 17. Fatigue growth onset for ENF specimen. 


41 


mm/cycle 



Figure 18. Delamination growth rate (Paris Law) for ENF specimen. 
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Figure 19. Delamination growth rate behavior for ENF specimen for applied cyclic loading. 
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Figure 20. Delamination growth behavior for ENF specimen for applied cyclic loading. 
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Figure 21. Delamination growth rate behavior for ENF specimen for applied center deflection. 
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Figure 22. Stable delamination growth behavior for ENF specimen for applied center deflection. 
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Figure 23. Total delamination growth behavior for ENF specimen for applied cyclic loading. 
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Figure 24. Total delamination growth behavior for ENF specimen for applied center deflection. 

2200 


load 
Q, N 



1.0 1.5 2.0 

displacement w, mm 


Figure 25. Load-displacement behavior for an ENF specimen obtained from different models. 
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Figure 26. Computed strain energy release rate distribution across the width of an ENF specimen. 
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Figure 27. Computed critical load-displacement behavior for ENF specimen obtained 
from two-dimensional planar models (CPE4I) with different release tolerance settings. 
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Figure 28. Computed critical load-displacement behavior for ENF specimen obtained 
from two-dimensional planar models (CPE4I) with added contact stabilization. 
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Figure 29. Computed critical load-displacement behavior for ENF specimen obtained 
from two-dimensional planar models with added viscous regularization. 
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Figure 30. Critical displacement-propagation behavior for an ENF specimen. 
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Figure 3 1 . Computed critical load-displacement behavior for ENF specimen 
obtained from results using three-dimensional solid finite element models. 
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Figure 32. Computed critical load-displacement behavior for ENF specimen 
obtained from results using continuum shell finite element models. 
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Figure 33. Computed critical displacement-propagation for ENF specimen 
obtained from results using three-dimensional finite element models. 
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Figure 34. Computed critical load-displacement behavior for ENF specimen obtained 
from two-dimensional planar models with added contact stabilization. 
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Figure 35. Computed critical load-displacement behavior for ENF specimen obtained 
from two-dimensional planar models with added viscous regularization. 
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Figure 36. Computed critical load-propagation for ENF specimen obtained from 
results using two-dimensional planar finite element models for applied center deflection. 
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Figure 37. Computed critical load-displacement behavior for an ENF specimen 
obtained from three-dimensional models subjected to an applied center load Q. 
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Figure 38. Computed critical load-propagation for ENF specimen obtained from 
results using three-dimensional solid finite element models subjected to an applied center load, Q. 
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Figure 39. Fracture surface of typical ENF fatigue specimen /(57- 
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a. Initial delamination front shape (Bottom surface of FE model in Figure 6). 
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Figure 40. Computed delamination front shape for ENF specimen (FE model in Figure 6). 


53 





a. Initial delamination front shape ( Bottom surface ofFE model). 



b. Detail of delamination initiation at c. Detail of delamination propagation, 

edges. 

Figure 41. Computed delamination front shape for ENF specimen (FE model with Ac/ =0.5 mm). 
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Figure 42. Computed delamination front shape for ENF specimen (FE model in Figure 7a). 


54 









time t, s 


Figure 43. Computed cyclic fatigue loading for ENF specimen. 
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Figure 44. Computed energy release rate for ENF specimen. 
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Figure 46. Computed delamination onset and growth obtained for different element length, Aa. 
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Figure 47. Computed delamination onset and growth for different applied cyclic displacements. 
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Figure 48. Computed delamination onset and growth obtained from full three-dimensional models 

(Aa=1.0 mm, Fig. 6). 
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Figure 49. Computed delamination onset and growth obtained from full three-dimensional models. 



Figure 50. Computed delamination onset and growth obtained for different values of t . 
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Figure 51. Computed delamination onset and growth obtained for different element lengths, Aa. 
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Figure 52. Computed delamination onset and growth for different applied cyclic loads. 
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Figure 53. Computed delamination onset and growth obtained from 
full three-dimensional models (Aa=1.0 mm, Fig. 6). 
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Figure 54. Computed delamination onset and growth obtained from 
full three-dimensional models ( Aa=0.5 mm). 
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b. Detail of delamination initiation at edges, c. Detail of delamination growth. 



d. Detail of delamination growth. e. Detail of final front shape. 

Figure 55. Computed delamination front shape for ENF specimen with one Fourier term 

(FE model in Figure 6). 
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b. Detail of delamination initiation at edges. c. Detail of delamination growth. 



d. Detail of delamination growth. e. Detail of final front shape. 


Figure 56. Computed delamination front shape for ENF specimen with 50 Fourier terms 

(FE model in Figure 6). 
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a. Detail of delamination growth. b. Detail of final front shape. 

Figure 57. Computed delamination front shape for ENF specimen with one Fourier term 

(Aa=0.5 mm). 



a. Detail of delamination growth. b. Detail of final front shape. 

Figure 58. Computed delamination front shape for ENF specimen with 50 Fourier terms 

(Aa=0.5 mm). 




a. Detail of delamination growth. b. Detail of final front shape. 

Figure 59. Computed delamination front shape for ENF specimen (A a=0.5 mm, model Fig. 7a). 
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delaminated section 


a. Initial delamination front shape ( Bottom surface ofFE model Aa=l .0 mm). 



d. Detail of delamination growth. e. Detail of final front shape. 


Figure 60. Computed delamination front shape for ENF specimen (FE model in Figure 6). 


64 




a. Detail of delamination growth. 


b. Detail of final front shape. 


Figure 61. Computed delamination front shape for ENF specimen with one Fourier term 

(Aa=0.5 mm). 




a. Detail of delamination growth. b. Detail of final front shape. 

Figure 62. Computed delamination front shape for ENF specimen with 50 Fourier terms 

(A a=0.5 mm). 




a. Detail of delamination intitiation at edges, b. Detail of final front shape. 


Figure 63. Computed delamination front shape for ENF specimen (A a=0.5 mm, model Fig. 7a). 
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Figure Al. Deformed 2D finite element model of an ENF specimen. 
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Figure A2. Mixed mode fracture criterion. 
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Figure A3. Delamination growth onset for mode II for IM7/8552 [6], 
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Figure A4. Delamination growth rate ( Paris Law) for IM7/8552. 
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time t, s 


Figure A5. Amplitude curve ( relative load magnitude ). 
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Figure A6. Cyclic fatigue loading input. 
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